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The dynamics of a nonconservative Gross-Pitaevskii equation for trapped atomic systems with at- 
tractive two-body interaction is numericaUy investigated, considering wide variations of the noncon- 
servative parameters, related to atomic feeding and dissipation. We study the possible limitations 
of the mean field description for an atomic condensate with attractive two-body interaction, by 
defining the parameter regions where stable or unstable formation can be found. The present study 
is useful and timely considering the possibility of large variations of attractive two-body scattering 
lengths, which may be feasible in recent experiments. 
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I. INTRODUCTION 

The stability of the condensed state is governed by the 
nature of the effective atom-atom interaction: the two- 
body pseudopotential is repulsive for a positive s— wave 
atom-atom scattering length and it is attractive for a 
negative scattering length [1]. The ultra-cold trapped 
atoms with repulsive two-body interaction undergoes a 
phase-transition to a stable condensed state, in several 
cases found experimentally, as for *^Rb [2], ^^Na [3] and 
[4]. However, a condensed state of atoms with nega- 
tive s— wave atom-atom scattering length (as in case of 
^Li [5]) would be unstable, unless the number of atoms N 
is small enough such that the stabilizing force provided 
by the zero-point motion and the harmonic trap over- 
comes the attractive interaction, as found on theoretical 
grounds [6,7]. Particularly, in the case of ^Li gas [5], for 
which the s— wave scattering length is a = —14.5 ± 0.4 
A, it was experimentally observed that the number of al- 
lowed atoms in the Bose condensed state was limited to 
a maximum value between 650 and 1300, a result consis- 
tent with the mean-field prediction [6], where the term 
proportional to the two-body scattering length (negative) 
dominates the nonlinear part of the interaction. 

More recently, the maximum critical number of atoms 
for Bose-Einstein condensates with two-body attractive 
interactions have been deeply investigated by the JILA 
group, considering experiments with *^Rb [8]. They have 
considered a wide tunning of the scattering length, a, 
from negative to positive, by means of Feshbach reso- 
nance [9,10], and observed that the system collapses for 
a number of atoms smaller than the theoretically pre- 
dicted number. Their experimental results, when com- 
pared with theoretical predictions for spherical traps, 
show a deviation of up to 20% in the critical number. 
More precisely, it was shown in Ref. [11] that part of this 
discrepancy is due to the non-spherical symmetry that 
was considered in Ref. [8]. Such a deviation can also be 
an indication of higher order non-linear effects that one 
should take into account into the mean-field description. 
In Ref. [12], it was considered the possibility of a real 



and positive quintic term, due to three-body effects, in 
the Gross-Pitaevskii formalism. A negative quintic term 
would favor the collapse of the system for a smaller crit- 
ical number of atoms, as verified in the JILA's experi- 
ments. However, the real significance of a quintic term 
in the formalism is still an open question. 

Our main motivation in the present work is to ana- 
lyze the dynamics represented by an extension of the 
mean-field or Gross-Pitaevskii approximation, with non- 
conservative imaginary terms that are added to the real 
part of the effective interaction, the two-body nonlinear 
term with a spherically symmetric harmonic trap. For 
the imaginary part, the interaction is a combination of a 
linear term, related to atomic feeding, and a quintic term, 
due to three-body recombination, that is responsible for 
the atomic dissipation. This is an approximation that is 
commonly used to study the properties of Bose-Einstein 
condensed systems. We consider a wide variation of the 
nonconservative parameters, in particular motivated by 
the actual realistic scenario, that already exists, of alter- 
ing experimentally the two-body scattering length [10]. 
As it will be clear in the following, this possibility will 
lead effectively to a modification of the dissipation pa- 
rameter. By changing the absolute value of the scattering 
length, from zero to very large absolute values, one can 
change in an essential way the behavior of the mean-field 
description. As it will be shown from the present numer- 
ical approach, the results for the dynamical observables 
of the system can be very stable (solitonic-type) or very 
unstable (chaotic- type); the characteristic of the results 
will depend essentially on the ratio between the noncon- 
servative parameters related to the atomic feeding and 
dissipation. 

In the next section, we review the formalism. The 
main results are presented in section III, followed by our 
conclusions in section IV. 
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II. MEAN FIELD APPROXIMATION 

The mean field approximation has shown to be appro- 
priate to describe atomic Bosc-Einstcin condensation of 
a dilute gas of atoms confined by a magnetic trap [13]. 
In the case of positive scattering length, we have a very 
good agreement with experimental data, as the thermal 
cloud is practically absent (removed by cooling evapora- 
tion) and almost all the particles are in the condensed 
state. In this case, the mean field approximation results 
in a nonlinear Schrodinger equation (NLSE) known as 
Gross-Pitaevskii equation (GPE) . If wc have N particles 
trapped in a spherical harmonic potential, this equation 
is given by 



where is the dissipation parameter, due to three-body 
collisions, and G-y is a parameter related to the feeding of 
the condensate from the thermal cloud. The Eq. (2) was 
first suggested in Ref. [14] to simulate the condensation 
of ^Li. 

In order to recognize easily the physical scales in 
Eqs. (1) and (2), it is convenient to work with dimen- 
sionless units. By making the transformations 
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^r *, (1) we obtain the radial dimensionless s— wave equation: 



where 4* = \l/(r, f), the wave-function of the condensate, 
is normalized to the number TV, m is the mass of a single 
atom, uj is the angular frequency of the trap, and a is the 
two-body scattering length. 

In this work, wc have concentrated our study in the 
interesting dynamics that occurs when the scattering 
length is negative (a = — |a|). In this case, it is well 
known that the system is unstable without the harmonic 
trap; and the trapped system has a critical limit Nc in 
the number of condensed atoms. The mean- field approx- 
imation has also shown to be reliable in determining the 
critical number of particles and even collapse cycles in 
the condensate [5,14,15]. Actually, systems with attrac- 
tive two-body interaction are being intensively investi- 
gated experimentally [8] , by using the so-called Feshbach 
resonance [9,10]. The scattering length can be tuned over 
a large range by adjusting an external magnetic field (for 
more details, see Ref. [16]). Here, we are interested in the 
dynamics of a realistic system, where wc add two non- 
conservative terms: one (linear) related to the atomic 
feeding from the non-equilibrium thermal cloud; and an- 
other, dissipative due to three-body recombination pro- 
cesses (quintic). It is true that other dissipative terms can 
also be relevant for an arbitrary trapped atomic system, 
as a cubic one, that can be related with dipolar relax- 
ation or with an imaginary part of the two-body scat- 
tering length. However, in order to simplify the study 
and better analyze the results, we restrict our considera- 
tions to the case that we have just one parameter related 
with the feeding and another related with dissipation. 
We have considered only the three-body recombination 
parameter for dissipation also motivated by the observa- 
tion that, for higher densities, this term dominates the 
two-body loss [17]. So, for the generalization of Eq. (1), 
we add the imaginary terms in the interaction, such that 
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As ^(r, t) is normalized to the number of atoms N{t) in 
Eq. (2), the corresponding time-dependent normalization 
of '^{x, t) is given by the reduced number n(r): 



dx\^{x,T)\'^ = n{T) =2N{t)\a\ 
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The nonconservative GPE (6) is valid in the mean-field 
approximation of the quantum many-body problem of 
a dilute gas, when the average inter-particle distances 
are much larger than the absolute value of the scattering 
length; and also when the wave-lengths arc much larger 
than the average inter-particle distance. The nonconser- 
vative terms are important when the condensate oscil- 
lates, fed by the thermal cloud, while losing atoms due 
to three-body inelastic collisions, which happen mainly 
in the high density regions. 

In order to verify the stability and the time evolution 
of the condensate, as observed in Refs. [18], two possible 
relevant observables are the number of particles normal- 
ized by the critical number of atoms of the static case 
{N{t)/Nc) and the mean square radius (msr). 
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where X = \/{xHt)) is the dimensionless root mean 
square radius (rmsr). In our analysis of stability, we cal- 
culate the time evolution of these quantities. We explore 
several combinations of the dimensionless nonconserva- 
tive parameters £, and 7. We first consider the case in 
which the atomic feeding is absent or when its param- 
eter is smaller than the atomic dissipation parameter. 
Next, we explore variations of both parameters of about 
five orders of magnitude, from 10~^ to 10~^. This wide 



2 



spectrum includes the parameters considered by Kagan 
et al. [14], as well as other combinations that can be 
considered more realistic due to recent experimental re- 
sults [19]. 

Actually, the relevance of a wider relative variation of 
the nonconservative parameters 7 and ^, presented in 
Eq. (6), can be better appreciated in face of the experi- 
mental possibilities that exist to alter the two-body scat- 
tering length [10]. As one should note from Eq. (4), any 
variation of the scattering length will also affect the effec- 
tive dissipation parameter ^ and, consequently, its rela- 
tion with the feeding parameter 7. This implies that, by 
changing the value of the scattering length, from positive 
to negative, and from zero to very large absolute values, 
one can change in an essential way the behavior of the 
mean-field description. In the present work, we are con- 
cerned with negative two-body scattering length, where 
the collapsing behavior of the Eq. (6) shows a very inter- 
esting dynamical structure. Even considering the possi- 
ble limitations on the validity of the mean-field approach 
after the first collapse (in cases of parametrization where 
it can occur), it is worthwhile to verify experimentally 
the behavior of a system in such a situation, by vary- 
ing |a|. At least, one can verify how far the theoretical 
description can be qualitatively acceptable. 

As already verified for systems with attractive interac- 
tion, as the ''Li, it has been possible, via the mean-field 
approach, to describe properties like the critical number 
of atoms in the condensate and growth and collapse cy- 
cles [5,14,15]; besides, in the long time evolution, for cer- 
tain sets of parameters, the calculations have also shown 
the presence of strong instabilities of the condensate, 
with signals of spatio-temporal chaotic behavior. 

In order to characterize a chaotic behavior, it is neces- 
sary to show that the largest Lyapunov exponent related 
with the solutions of the equation is positive. We follow 
the criterion used by Deissler and Kaneko [20] to char- 
acterize spatio-temporal chaos. This criterion prescribes 
that the largest Lyapunov exponent for the system, in 
an arbitrary time interval, is obtained by plotting the 
logarithm of a function ^, which is defined by 

/ foo \ 1/2 

C(r)^^ mx,T)\^dx^ . (9) 

5^{x,t) will give us the separation between two nearby 
trajectories; it is obtained in the following form: we 
numerically evolve in time an initial $0(2;), obtaining 
$(x,t). Independently, we evolve $0(2;) + and get 
$'(x,r), where e(x) is a very small random perturba- 
tion, r) is given by $'(a;, r) — $(a;, r). The chaotic 
behavior is characterized by a positive slope of ln(^(T), 
which gives the largest Lyapunov exponent [20]. 



III. NUMERICAL RESULTS 

In the next, we present the most significant results that 
characterize the time evolution of the normalized number 
of particles ( N{t)/Nc ), the dimensionless mean square 
radius (a;^), and, in order to characterize the stability of 
the system, the function related to the largest Lyapunov 
exponent. Further, we present a representative case of 
the phase-space for the root mean square radius. We 
have studied a wide region of parameters 7 and ^, cover- 
ing about five orders of magnitude, from 10~^ to 10~^, 
including the case with no feeding (7 = 0). 

In order to have a clear and useful map of the regions 
where one should expect stable results, as well as regions 
with instabilities or chaos, we summarize the present nu- 
merical results in a diagrammatic picture that relates 
these two nonconservative parameters. In general, it is 
expected that the system is more stable when the param- 
eter related to the feeding of atoms (7) from the thermal 
cloud is significantly smaller than the parameter related 
to the dissipation (^). However, it is interesting to find 
out the region of parameters where this transition (from 
stable to unstable results) occurs. Analysis of experimen- 
tal results can provide a test to the present mean-field 
description in case of negative two-body interaction. As 
previously observed, we are considering dimensionless ob- 
servables and parameters. For any realistic comparison 
with experimental parameters, one should convert 7 and 
^ to the parameters and G^, as given in Eq. (4). 

The numerical solutions of Eq. (6) were obtained by 
applying the semi-implicit Crank-Nicolson algorithm for 
nonlinear problems, as implemented in Ref. [18]. This 
method is stable and, therefore, very convenient and re- 
liable to treat time-dependent non-linear partial differ- 
ential equations. The initial condition for the number of 
atoms A'' in the condensate was such that N{0)/Nc = 
n(0)/r?,c = 0.75. The evolution of the observables have 
been extended up to r = = 500. 

In general, as expected, the smaller is the dissipation 
parameter, the longer is the life of the condensate. The 
mean sqiiare radius presents an oscillatory behavior while 
one increases ^. One observes that, in the regime of small 
feeding (7 < 10"''), the extended Lyapunov presents no 
positive slope. For larger values of 7, from ~ 10~^ and 
10^^, we have studied a few cases where the interplay 
between the nonconservative behaviors are significant. 

In Fig. 1, we show the dynamical behavior of the num- 
ber of atoms for 7 = 10^^ and several values of ^; and, in 
the Fig. 2, the corresponding time evolution of {x{t)'^). 
We realize an interesting behavior, that occurs when the 
dissipation is larger than the feeding process: there are 
solutions of stability or dynamical equilibrium between 
both nonconservative processes. This phenomenon was 
already discussed in Ref. [21], for a few values of the dis- 
sipation and feeding parameters, using the time depen- 
dent variational approach and also the Crank-Nicolson 
method. In the present work, we observe a wide region of 
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parameters where it is possible the formation of autosoli- 
tons [21]. However, when the feeding process is much 
larger than the dissipation, of about one or more orders 
of magnitude, we can also observe chaotic behaviors. See, 
for example, the case with ^ = 10"'^. 

The time evolution of the number of particles, repre- 
sented in Fig. 1, shows a collapse for r «30, followed by 
several other collapses, with the number of particles go- 
ing above the critical limit Nc- So, after a sequence of 
collapses, the critical limit for the number of particles is 
no more followed, as already shown in Ref. [18] . 
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FIG. 1. Time evolution of the number of condensed atoms 
N, relative to the critical number Nc, for a set of values of 
the dissipativc parameter ^ (as shown inside the frame) , with 
the feeding parameter 7 = 10~^. All the quantities are in 
dimensionless units, as given in Eqs. (3)- (4). 

The corresponding time evolution of {x(t)'^) is shown 
in the upper frame of Fig. 2. We observe that, follow- 
ing each collapse, after the shrinking of the system, the 
radius is multiplied by a large factor, with indication of 
being populated by radial excited states. In the lower 
frame of Fig. 2, we can observe the corresponding tran- 
sition from the stable region (where the system finds the 
equilibrium at a fixed value of the radius, correspond- 
ing to autosoliton formation) to the unstable region. As 
shown, the instability starts to occur when ^ = 2 x 10~^, 
and it can be developed to a spatio-temporal chaos. The 
chaotic behavior can be verified through the Deissler and 
Kaneko criterion [20]. 
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FIG. 2. Time evolution of the dimensionless mean-square- 
-radius, (a;(r)^), for the feeding parameter 7 = 10~^. The 

results arc given for a set of values of the dissipativc parameter 
^, in the lower frame (shown inside). A specific case, for 



10" 



(much smaller than 7), is isolated in the upper 



frame, where one can observe the behavior of {x{t)'^) after 
the collapse. All the quantities are in dimensionless units, as 
given in Eqs. (3)-(4). 

In Fig. 3 we illustrate the application of the Deissler 

and Kaneko criterion to the system given by Eq. (6), for 
a fixed value of the feeding parameter 7 = 0.01, and 
a set of values of the dissipation parameter ^. It was 
plotted the time evolution of the function ln(C), where 
C is given by Eq. (9), following the prescription given 
in Ref. [20] to obtain the largest Lyapunov exponents 
for the system. Within this prescription, the system be- 
comes chaotic when In(^) has a positive slope. As shown 
in Fig. 3, this clearly occurs, for example, when ^ = 10~^. 
In case of ^ = 10~^ we note a much faster increasing in 
In(^), with an observed saturation that happens due to 
the fact that such function has reached the maximum sep- 
aration between the trajectories. The saturation prop- 
erties is also verified when studying chaotic behaviors in 
ordinary differential equations [22]. The plot of ln(C) cor- 
responds to the same value of 7 (= 10~^) used in Figs. 1 
and 2. As shown, a clear characterization of chaotic be- 
haviors starts to occur only for values of the dissipation 
parameter ^ much smaller than 7. In the cases presented 
in Fig. 3, for ^ < 10"^. 
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FIG. 3. Time evolution of Ini^, related to the separation 
between two nearby trajectories [See Eq. (9)], for 7 = 10~^ 
and a set of values of ^ indicated inside the figure. All the 
quantities are in dimensionless units, as given in Eqs. (3)- (4). 

In Fig. 4, wc present another significative illustration 
of chaotic behavior, through the phase-space behavior 
of the mean-square-radius, considering one case that was 
characterized as chaotic by using the Deissler-Kaneko cri- 
terion. We have plotted in this figure the root mean- 
square-radius phase space for the case vi^ith 7 = 0.01 and 
^ = 10~^. The irregular behavior of the trajectories, 
observed in Fig. 4, with the classical strange attractors 
being observed, clearly resembles chaos. This behavior 
is similar to the chaotic behavior observed in ordinary 
cases [22]. 
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FIG. 4. Phase-space for the root mean-square-radius, in di- 
mensionless units [dX{T)/dt versus X(r)], considering a col- 
lapsing case that leads to chaos. The dimensionless noncon- 
servative parameters are ^ = 10"'' and 7 = 10~^, and the 
time evolution was taken up to t = r/w = SOO/w. 



sented results, wc should note that, in order to observe 
unstable chaotic behaviors, the dissipation must be much 
smaller than the feeding parameter. 

In a diagrammatic picture, given in Fig. 5, we resume 
our results. We show the relation between the two non- 
conservative parameters, ^ and 7, in order to characterize 
the parametric regions where one should expect stabil- 
ity or instability in the solutions for the Eq. (6). The 
stable results of the Eq. (6) arc represented by bullets; 
the nonstable results that clearly present positive slope 
for InC(T) (chaotic behavior) are represented by empty 
squares; with x , we show other intermediate nonstable 
results, in which the characterization of chaotic behav- 
ior was not so clear, through the Deissler-Kaneko crite- 
rion. In this figure, in order to observe the approximate 
consistency of the numerical results, we also include the 
variational analysis presented in Fig.l of Ref. [21], rep- 
resented by the dashed-line. It is separating the stable 
region (upper part) from the unstable one (lower part). 
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As a general remark that one can make from the pre- 



FIG. 5. Diagram for stability, according to the criterion of 
Ref. [20] given by Eq. (9), with results for the equation (6), 
considering the dimensionless non-conservative parameters 7 
and ^. Between the unstable results, represented with x and 
squares, the chaotic ones are identified with squares. The sta- 
ble results arc represented by bullets. Two dotted guide-lines 
are splitting the regions. The dashed line split the graph in 
two regions according to a variational approach (see Ref. [21]); 
in the upper part the results are stable, in the lower, unstable. 

We should note that, in section V of Ref. [23], it was 
also considered the dynamics of growth and collapse, with 
non-conservative terms related to feeding (70) and dissi- 
pation (71 and 72) in a specific example. For the dissipa- 
tion they have also considered a term related to dipolar 
relaxation, given by 72. Here, in our systematic study of 
the regions of instability, wc took into account previous 
experimental [8] observations that the dominant process 
for the dissipation is the three-body recombination. By 
comparing the parameters of Ref. [23] with the parame- 
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ters that we have used, and observing that our parameter 
^ should be related to both dissipation parameters used 
in Ref. [23] (7 = 70 = 2.6 x 10"^^ ^ _ 
verify from the results given in Fig. 5 that the model 
of Ref. [23] is inside the intermediate region, where the 
system is unstable, without a clear signature of chaos. 

IV. CONCLUSIONS 

In summary, we have studied the dynamics associ- 
ated with the extended nonconscrvativc Gross-Pitaevskii 
equation for a wide region of the dimensionless noncon- 
servative parameters, ^ and 7, that, respectively, are 
related to atomic dissipation and feeding in a trapped 
atomic condensed system. We consider systems with at- 
tractive two-body interaction in a spherically symmet- 
ric harmonic trap. In Fig. 5, we resume our results, by 
mapping the space of 7 versus ^, showing the regions 
of equilibrium and the regions of instability, as well as 
the regions where we are able to characterize chaotic be- 
haviors, using a criterion given in Ref. [20]. It was also 
confirmed that chaotic behaviors occur mainly when 7 is 
big enough and 7/^ is large (at least, when 7 is one or 
two orders of magnitude larger than ^). A wide variation 
of the nonconscrvativc parameters was analyzed, in par- 
ticular motivated by the actual realistic scenario, that 
already exists, of altering experimentally the two-body 
scattering length [10]. By changing the absolute value of 
the scattering length, one can change in an essential way 
the behavior of the mean-field description. 
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